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ABSTRACT 

^  Recently  the  authors  showed  the  convergence  of  a  class  of  vortex  methods 
for  incompressible ,  inviscid  flow  in  two  or  three  space  dimensions.  These 
methods  are  based  on  the  fact  that  the  velocity  can  be  determined  from  the 
vorticity  by  a  singular  integral.  The  accuracy  of  the  method  depends  on 
replacing  the  integral  kernel  with  a  smooth  approximation.  The  purpose  of 
this  note  is  to  construct  smooth  kernels  of  arbitrary  order  of  accuracy  which 
are  given  by  simple,  explicit  formulas. 
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SIGNIFICANCE  AMD  EXPLANATION 


Vortax  methods  simulate  incompressible,  inviacid  flow  by  a  system  of 
ordinary  differential  equations  for  the  paths  of  representative  particles  in 
the  fluid.  They  have  the  advantage  that  the  computational  elements  are  auto¬ 
matically  concentrated  in  the  region  of  vorticity  and  errors  like  the 
numerical  diffusion  of  difference  methods  are  avoided.  The  use  of  modified 
velocity  kernels  ensures  the  accuracy  and  stability  of  the  method.  Explicit 
formulas  for  these  kernels  make  it  possible  to  implement  this  method  as 
directly  and  efficiently  aa  if  no  smoothing  were  used. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


EXPLICIT  SMOOTH  VELOCITY  KERNELS  FOR  VORTEX  METHODS 

*  ** 

J.  T.  Beale  and  A.  J.  Majda 

Recently  the  authors  proved  the  convergence  of  a  class  of  vortex  methods 
for  the  simulation  of  incompressible/  inviscid  flow  in  two  or  three  space 
dimensions  without  boundaries  (see  [1,2]).  The  principle  of  such  methods  is 
to  reduce  the  calculation  of  the  flow  to  a  system  of  ordinary  differential 
equations  for  the  paths  of  representative  particles.  The  velocity  field  is 
given  by  the  convolution  of  a  singular  kernel  with  the  vorticity 
distribution.  In  the  class  of  methods  treated,  the  stability  and 
discretization  error  are  controlled  by  a  distortion  of  the  nearby  interaction 
of  the  particles.  This  is  accomplished  by  convolving  the  integral  kernel  with 
a  smooth  approximation  to  the  delta  function.  The  choice  of  this  function 
determines  the  order  of  accuracy  of  the  method  for  smooth  flows.  The  purpose 
of  this  note  is  to  produce  a  class  of  functions  which  lead  to  smooth  kernels 
given  by  simple,  explicit  formulas.  With  these  choices,  the  method  can  be 
implemented  with  essentially  no  more  effort  than  would  be  necessary  if  the 
original  kernel  were  used.  In  each  case  treated  here,  kernels  of  high  order 
accuracy  are  obtained  easily  from  ones  of  low  order  by  scaling  or  other 
modification. 

A  general  description  of  various  vortex  methods  in  use  can  be  found  in 
the  survey  article  of  Leonard  [6].  A  summary  of  theoretical  results  is 
contained  in  [7],  and  complete  proofs  are  given  in  [1,2,5).  Numerical 
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experiments  of  Held  and  Del  Prete  [4]  for  two-dimensional  methods  as  in  [2] 
are  in  general  agreement  with  the  predictions  of  the  convergence  results,  and 
further  experiments,  including  tests  of  higher  order  accuracy,  are  currently 
under  way.  A  three-dimensional  method  similar  to  that  of  [1],  but  with 
significant  differences,  has  been  used  by  Chorin  [3]. 

Vortex  methods  are  based  on  the  fact  that,  for  Incompressible  flow,  the 
velocity  is  determined  from  the  vorticity  by  a  convolution, 

u(z,t)  “  (K*o>) <s,t)  ■  /  X(z-x' ) uK*' ,t)dz*  .  (1) 

(We  will  use  notation  consistent  with  (1,2].  This  formula  is  interpreted 
differently  in  two  or  three  dimensions)  see  below.)  In  the  methods  of  [1,2], 
as  in  earlier  work,  the  velocity  kernel  K  is  replaced  by 

K5  -  K  *  ♦g,  ♦gts)  -  6"%<z/6)  .  (2) 

Here  M  *  2  or  3  is  the  space  dimension  and  S  isa  parameter  to  be  chosen 
in  conjunction  with  the  linear  spacing  h  of  the  particles  introduced  at  time 
zero.  The  smoothing  of  the  kernel  by  the  function  ^  can  be  interpreted  as 
the  approximation  of  the  vorticity  distribution  by  a  sum  of  "blobs”  of 
prescribed  shape  (see  [5]  or  [7]). 

We  will  choose  the  function  1>  subject  to  the  conditions 

(i)  ^  is  smooth  and  rapidly  decreasing,  i.e., 

|dS(*)I  <  Cgj(1  +  |z|Vj 

for  every  multi-index  8  and  every  integer  j) 

(ii)  /  ♦(z)dN*  -  1  i 

(iii)  /  z%(*)dNz  -  0,  1  <  181  4  p-1  , 

where  p  is  an  integer. 

The  results  of  [1,2]  imply  that  vortex  methods  satisfying  (i)  -  (iii) 
converge  provided  the  relation  between  6  and  h  is  properly  chosen.  If 
4  -  hq,  q  any  fixed  number  with  0  <  q  <  1,  the  error  is  of  the  order  of 


6P  *  h^*3,  i.e.  the  method  is  essentially  pth  order.  Our  object  here  is  to 
choose  so  that  Kg  has  a  simple  expression  consistent  with  these 

requirements.  As  we  shall  see,  choices  of  with  p  *  2  easily  lead  to 
choices  with  p  >  4. 

Condition  (i)  implies  that  the  Fourier  transform  of  as  well  as  ^ 

itself,  is  smooth  and  rapidly  decreasing.  We  will  always  take  *  4*(r),  r  * 

| z | .  In  this  case  (iii)  holds  by  symmetry  for  | 6|  odd,  so  that  p  may  be 
assumed  even.  For  radial  f,  (iii)  is  always  satisfied  with  p  *  2.  This 
set  of  conditions  is  more  stringent  than  those  in  the  general  theory  of  [1,2]. 
(In  the  earlier  language  we  are  assuming  ^  to  be  in  the  class  Fe  S  . ) 

The  condition  (i)  can  be  relaxed  somewhat  to  allow  a  which  is  not  very 
smooth  at  z  *  0.  Indeed,  our  simplest  choice  in  three  dimensions  has  this 
property.  With  a  weaker  condition  replacing  (i),  a  similar  convergence  result 
holds,  but  q  is  restricted  to  an  interval  0  <  q  <  q0  for  some  q0  <  1. 

(See  [1]  for  precise  statements.) 

Two-dimensional  Flows 

In  the  2-D  case  the  vorticity  is  the  scalar  function  <■>  =*  u  -  u 

* ,  x  ■ ,  y 

The  distinguishing  property  of  two-dimensional  flows  is  that  the  vorticity  is 
conserved  along  particle  paths: 

»  +  u  •  V»  *  0  .  (3) 

Suppose  an  initial  velocity  field  is  prescribed  with  vorticity  nonzero 

only  within  a  bounded  set.  To  simulate  the  flow,  we  cover  this  set  with  a 
square  grid  of  size  h  and  introduce  a  particle  at  the  center  of  each 
square.  We  take  the  coordinates  of  a  typical  particle  to  be  ih,  where  i  « 

( i ^ ,  i 2 )  is  a  pair  of  integers;  the  ith  particle  is  assigned  the  vorticity 
■  wQ(ih).  To  compute  approximate  paths  of  the  particles,  we  discretize 
(1),  with  K  replaced  by  K.,  remembering  (3),  and  arrive  at  the  system  of 


ordinary  differential  equation 

f  ^  2 

z.  -  )  K .  (z  -  z.)a»,h 
i  j  o  i  3  j 


zi(0)  *  ih  . 


The  area-preserving  property  of  incompressible  flow  is  used  implicitly  here. 
Once  the  z^s  have  been  determined,  an  expression  for  the  veloc 
be  obtained  by  setting 

u*\z,t)  -  l  Kfi(z-z .. )oyi2  . 

To  apply  this  method  it  is  best  to  have  an  explicit  formula  for  K_. 

6 

If  G  is  the  Green's  function  for  -V2,  G(z)  *  -(2»)  1  log  r,  then  with  z  * 


(x*y) 


K(z)  -  (8  ,  -3  )G  «  . 

y  X  „  2 

2  2*r 


( 2 )  -r 

A  natural  choice  of  +  is  the  Gaussian  ^  (r)  *  e  /w.  The  necessary 

conditions  (i)  -  (iii)  are  satisfied  with  p  *  2.  If  K.  **  K  *  i|>.  then 

0  0 

K.  -  (3  ,  -3  )G.  ,  G.  -  G  *  . 

o  y  x  o  6  T6 

Since  i|t  is  radial,  G^  is  also,  and 


V2G4  -  V2(G  *  ^6)  -  or 


;  V'W  -  -*«">  -  -  -7  ‘~eW 

*  6 


after  integration  we  have 


°r°S  ‘  2fr  (e'r2/{2  -  11  • 

The  constant  of  integration  is  determined  by  the  fact  that  G  must  be 

6 

smooth.  The  corresponding  velocity  kernel  may  now  be  found  from  (5)s 

k‘2,(.)  _  ,,  .  ,-tW,  . 

5  2  Hr2 


The  superscript  (2)  has  been  inserted  to  indicate  the  order  of  the  kernel. 

(4) 

Next  we  will  obtain  a  fourth  order  kernel  by  choosing  i1  m  i>  as  a 
combination  of  two  Gaussians  with  different  scalings, 

♦(4)<r)  -  Cl*<2>(r)  +  c2*(2)(r/a> 

where  a  is  arbitrary  except  that  a  +  1.  To  satisfy  condition  (ii)  we  must 


2 

c^  +  a  Cj  ■  1 

This  leaves  us  with  one  constraint  to  impose  moment  conditions.  Because  of 
symmetry,  condition  (iii)  will  hold  with  p  *  4  provided 

/  f/4^(r)r2  •  rdr  ■  0  . 


This  in  turn  holds  if 


C«  ^  3L  Cty  *  0  f 


(4)  (4) 

and  the  two  equations  determine  <fi  in  terms  of  a.  We  can  now  find  K ^ 

just  as  in  the  previous  case: 

2  .2  ?  2 .2 

_( 4).  .  (-v.x)  -r  /6  2  -r  /a  6  , 

K'  (z)  -  (1  -  c,e  -  c.a  e  )  . 

6  2irr2  1  2 

2 

For  example,  the  choice  a  "2  leads  to 


,  . ,  ,  ,  2  ,  ,2  2  ,2 

K<4)U)  .  (1  -  2e-r  /«  +  /2«  , 

6  2wr2 


-  (1  -  e"r  /26  HI  +  2e“r  /26  )  . 

2tr 


It  should  be  clear  that  hitler  order  kernels  can  be  constructed  by  adding 
further  terms  with  Gaussians  of  different  scalings  in  the  expression  for  4*. 
A  typical  sixth  order  kernel  is 


*«•>(.)  -  (i  - 1  *  2.-'2/2<s2  - 1  .-'2/“2,  . 

8  2wr2  3  3 


Of  course,  care  must  be  taken  in  evaluating  any  version  for  r  small,  since 

the  factor  due  to  the  smoothing  vanishes  at  z  *  0. 

Even  simpler  high  order  kernels  can  be  obtained  by  choosing  <(/  in  the 
_  2 

form  i>(r)  -  P(r)e”r  ,  where  P  is  a  polynomial  in  even  powers  of  r.  Then 
as  in  the  argument  above 

2  ,2 

rD  G, (r )  -  -«"2  /  rP{r/6)e"r  /<r<ir 
r  6  * 

2  2 

-  (2w)"1{Q(r/5)e“r  /S  -Q(0)}  , 

where  Q  is  another  even  polynomial  of  the  same  degree.  For  condition  (ii) 

to  hold  we  must  have  rD  G,(r)  -1/2*  as  r  ♦  so  that  Q(0)  -  1. 

r  o 


Finally 


and  according  to  (5), 


2  g2 

DrVz)  '  2«  (2(r/6)e"r  /r  “  ^  ' 

2.2 

K.(z)  -  K(z) ( 1  -  Q(r/6)e“r  '  )  . 


To  satisfy  the  moment  conditions  (iii)  we  need  to  have 


r2jP(r)e-r  r  dr  -  0, 
or  after  integrating  by  parts. 


1  <  j  <  (p-2 )/2  , 


r2j_1Q(r)e"r  dr  -  0  . 

The  moment  conditions  are  thus  reduced  to  linear  equations  for  the 
coefficients  of  Q.  With  p  *  4  we  find  Q(r)  -  1  -  r2,  corresponding  to 
the  fourth  order  kernel 


.(4) 


2  ..2 


(z)  -  {1  +  (-1  ♦  rW).-*  /6  }  . 


2*r 


For  order  p  a  polynomial  of  degree  p-2  is  sufficient.  For  example, 
with  p  ■  6  we  have 

Kj6)  (z)  -  K(z){l  -  Q(r/6)e"r 
with  -fi(r)  -  -1  +  2r  -  r*/2. 


-6- 


Three-dimensional  flow 

In  the  three-dimensional  eaaa  the  vorticity  u  -  7 Hi  is  a  vector 
quantity,  and  the  velocity  is  expressed  in  the  fora  (1)  by  the  Biot-Savart 
Law.  We  will  approximate  the  particle  paths  by  an  analogue  of  (4),  but  now 
the  vorticity  oust  be  updated  as  well  as  the  position.  A  natural  way  to  do 
this  is  to  use  the  direct  generalization  of  (3) 

u>t  +  u»Vu"<i»»7u  t 

this  is  essentially  the  method  used  by  Chorin  [3].  However,  the  derivative  on 
the  right  must  be  computed  by  differences  based  on  the  current  particle 
locations,  which  already  contain  errors.  To  avoid  this  difficulty,  we  use 
instead  a  Lagrangian  expression  for  the  vorticity  evolution  due  to  Cauchy. 

Let  a  be  the  Lagrangian  position  and  i  a  ♦  z  the  coordinate  mapping 
induced  by  the  flow.  Then 

u(z,t)  »  VtNa)  •  ftL(a) 

where  u)Q  is  the  initial  vorticity  and  z  ■  •*(<*).  Thus  the  vorticity  is 
carried  along  particle  paths  but  distorted  by  the  Jacobian  matrix  of  the 
flow.  Differentiating  in  t,  with  a  fixed,  we  have 

<*/t>  -  Vau(z,t)  •  «Q(o)  , 

* 

with  a  Lagrangian,  rather  than  Euler ian,  gradient  of  the  velocity.  This 

derivative  can  easily  be  implemented  by  a  difference  operator  on  the  initial 

grid.  We  are  thus  led  to  the  following  method!  given  current  values  of 

{z^},  {a>j},  we  can  compute 

uh(z,t)  ■  l  K.(z  -  z  .(t)  )S.(t)h3  . 

j  *  j  3 

The  updates  are  then  defined  by  the  system  of  ordinary  differential  equations 

M  Mr  % 

zj^  •  U  (zA,t>,  Z^tO)  " 

iz^t)  •  (i>0 ( ih ) ,  «1(0)  -  «0(ih>  , 
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r 


where  ?r‘  is  an  anti- symmetric  difference  operator  whose  order  is  at  least 
the  intended  order  of  accuracy.  (See  [1]  for  more  details.) 

The  three-dimensional  realisation  of  ( 1 )  is 

u  *  K  *  w  -  V*(G  *  u)  , 

where  G  is  the  Green's  function  for  -V2,  G  "  1/4 nr,  and  the  convolution  is 
componentwise.  As  before  we  set  Gg  -  G  *  ^  and  Kg  ■  K  *  4^.  It  is  easy 
to  see  that 

3G.  , 

(Kg  *  w)(s)  -  /  (If-s'l)  jfEfrj  *  uKs')ds' 

or  more  briefly 

V*'  -  aT  ftr  -  • 

Thus  Kg  will  have  a  simple  expression  provided  3Gg/3r  does  so. 

In  this  case  it  is  less  clear  how  to  choose  4  than  in  the  two- 
dimensional  case,  and  it  is  best  to  proceed  in  the  opposite  direction  from 
before.  For  simplicity  we  assume  at  first  that  6  -  1.  We  suppose  that 


3r 


f  (r) 
4wr2 


with  f  to  be  determined;  this  form  is  convenient  since  we  expect 
'*G1/3r  -  3G/3r  *  -1/4mr2  as  r  ♦  •.  Then 

-  V2G1  -  r"2Dr(r2DrG1} 


so  that 


♦(r) 


f '  (r ) 
4irr2 


(6) 


We  can  now  list  the  conditions  which  must  be  satisfied  by  f  so  that 
conditions  (i)  -  (iii)  hold  for  4*  Since  4  should  be  smooth  we  require 
(Cl)  f'(r)  is  a  smooth  function  of  r3; 

(C2)  f(r)  *  0(r3)  as  r  ♦  0. 

It  is  easily  seen  that  the  total  weight  of  4  is  f («),  and  our  last 


condition  is  therefore 


(C3 )  f(r)  ♦  1  as  r  +  *,  and  f'  is  rapidly  decreasing. 

If  (cl)  -  (c3)  hold,  then  defined  by  (6),  satisfies  (i)  -  (iii)  with 

p  »  2.  Two  choices  of  f  meeting  our  requirements  are 

__3  , 

f(r)  =  1  -  e  r  ,  f(r)“  tanh  rJ  , 

3  3 

corresponding  to  i|i(r)  =*  (3/4ir)e  r  and  i[i(r)  “  (3/4w)sech  r  ,  respectively. 
The  first  choice  is  simpler  and  is  analogous  to  the  Gaussian  function  in  two 
dimensions.  Actually,  in  this  case  (Cl)  does  not  hold  in  the  strictest  sense 
at  the  origin  because  f(r)  has  terms  r®  and  higher  in  the  Taylor 
expansion.  However,  the  general  theory  of  [1,2]  applies  to  this  choice,  and 
we  do  not  expect  the  difference  to  be  significant  in  practice.  Having 
chosen  f,  we  define  i|».  from  (2)  and  reverse  our  steps  to  find 


so  that 


^6  f  (r/6) 


3r 


4»r 


K.(z)  - r  f(r/5)s  x  . 

6  4xr3 

The  kernels  just  constructed  are  second  order  accurate  with  respect  to 
6.  If  f'  is  arbitrarily  smooth,  as  is  true  for  f(r)  *  tanh  r3,  then  we 
have  convergence  with  6  *  h^,  q  “  1-e,  and  the  errors  are  essentially  second 
order  in  h  as  well.  For  the  "cubic  Gaussian",  the  results  of  [1]  require 
q  <  5/8,  so  that  the  order  of  convergence  in  h  is  5/4  -  e.  (See  Theorem  1 
in  [1]  and  the  remarks  following  Theorem  2;  the  number  M  is  6  in  this 
case.  The  predicted  order  of  convergence  is  not  sharp  and  can  be  improved  at 
least  to  3/2  -  e.) 

To  obtain  kernels  with  p  *  4  we  can  combine  two  different  scalings  as 
(4) 

before.  Let  <|»  (r)  ■  c^r)  +  c2<|Kr/a),  where  ip  is  one  of  the  two 

(4) 

choices  specified  above.  The  conditions  that  $  have  weight  one  and 
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I 


satisfy  the  second  order  moment  conditions  lead  to  the  equations 


C1  +  »  c2  ■  1 


c1  +  a  c2  -  0 


Again  reversing  the  steps  we  have 


♦?’<*>  ■ 


+  c2a2f '  (r/a6)} 


3Gi 

3r 


(r) 


4»r  6 

_ 1 

4irr 


2  (c^tr/fi)  +  c2a  f(r/a6)} 


K6(z) 


— j  (c1f(r/6)  +  c2a3f(r/a4)}*x  . 


4rr 


For  example,  with  f(r)  »  1-e_r  or  tanh  r3  it  would  be  convenient  to 
-3 


choose  a 


2.  In  the  first  case  this  gives 


V*> 


4*r 


3  ..3  „  3  . .3 

t,  -r  /6  -2r  /6  , 

J  (1  -  -  c2e  }z  x 


so  that  only  one  exponentiation  is  necessary. 

An  alternative  method  can  be  used  to  produce  kernels  of  fourth  order  in 
6  from  the  ones  of  second  order  already  obtained.  Suppose  a  function  f(r) 
has  been  found  as  above  meeting  conditions  (Cl)  -  (C3).  He  will  choose 


f4(r)  «  c^tr)  +  c2rf’(r) 

with  appropriate  constants  and  check  that  the  corresponding  kernel 

1 


(7) 


4irr' 


—  f4(r/6)z* 


is  fourth  order. 

If  f  satisfies  (CD  and  (C2),  then  f4  does  also,  and  (C3)  will  hold 

provided  c^  ■  1.  We  need  to  impose  the  moment  condition  (iii)  with 

(4) 

|0|  *  2.  The  correspondence  between  f4  and  \|>  is  given  by  (6),  and  we 

•  (4 ) 

can  convert  the  condition  on  (1  to  a  similar  one  for  f4: 
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0  -  4w  Jj  *u,(r)r2  •  r2dr  -  /J  f^(r)r2dr 
”  £  gj(r)r2dr  -  -2  £  g4(r)r  dr  , 

where  g4(r)  ■  f4<r)  -  1,  so  that  g4<“)  ■  0.  If  we  now  substitute  from  (7) 
with  f  »  1  +  g  and  »  1,  our  condition  becomes 

0  -  rQ  {g<*>*  +  c2g*(r)r2}dr 
-  (1  -  2c2)  g(r)r  dr  , 

after  integrating  by  parts.  This  holds  if  c2  -  j,  and  therefore  the  choice 

f4(r)  -  f (r)  j  rf'(r) 

satisfies  all  our  requirements. 

we  can  apply  this  result  directly  to  the  two  choices  of  f  made  above. 
-r3 

If  f(r)  *  1  -  e  r  ,  then 

3 

f4(r)  »  1  ♦  (-1  +  |  r  )e”r  . 

For  f(r)  *  tanh  r3,  we  have 

f4(r)  «  tanh  r3  +  j  t3  sech2  r3 
or 

f4(r)  -  T  +  j  r3{1  -  T2),  T  »  tanh  r3  . 

Again,  either  of  these  two  methods  can  be  extended  to  produce  higher 
order  kernels.  Similar  arguments  could  be  used  in  the  2-0  case  and  would 
lead  to  the  Gaussian  kernels  found  before  and  additional  ones  as  well. 

Although  convergence  proofs  for  vortex  methods  have  been  given  only  in  the 
circumstances  of  [1,2,5],  the  smooth  kernels  obtained  here  could  be  applied  to 
other  versions  of  these  methods,  for  example  the  three-dimensional  method  of 
[3].  It  should  be  noted,  however,  that  for  a  vorticity  distribution  on  a  set 
of  lower  dimension,  such  as  an  interface  between  two  potential  flows,  the 
formulas  for  modified  kernels  are  somewhat  different.  Of  course,  optimal 
choices  of  the  parameters  will  have  to  be  determined  through  detailed 
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numerical  experiments 
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